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Abstract 

We have studied the hadronic interaction for the calculation of the atmospheric neutrino flux 
by summarizing the accurately measured atmospheric muon flux data and comparing with simu- 
lations. We find the atmospheric muon and neutrino fluxes respond to errors in the vr-production 
of the hadronic interaction similarly, and compare the atmospheric muon flux calculated using the 
HKKM04 [l| code with experimental measurements. The /i"*" + fj,~ data show good agreement 
in the 1~30 GeV/c range, but a large disagreement above 30 GeV/c. The ratio shows 

sizable differences at lower and higher momenta for opposite directions. As the disagreements are 
considered to be due to assumptions in the hadronic interaction model, we try to improve it phe- 
nomenologically based on the quark parton model. The improved interaction model reproduces the 
observed muon flux data well. The calculation of the atmospheric neutrino flux will be reported 
in the following paper 2|. 
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I. INTRODUCTION 



Evidence of neutrino oscillations was found in the atmospheric neutrino data observed 
with Super-Kamiokande . Atmospheric neutrinos are still a powerful tool to study neutrino 
oscillations, and the overall uncertainties in the observed data are becoming smaller ^. It is 
highly desirable to predict the absolute flux value and ratios among different kind of neutrinos 
precisely, and to understand their "systematic" uncertainties. Note, atmospheric neutrino 
experiments cover a wide L/E{= [neutrino flight length]/ [neutrino energy]) range, over 
four orders of magnitude 15|, which is much wider than accelerator neutrino beam experi- 

n 

ments such as K2K [6|. Atmospheric neutrino experiments are complementary to accelerator 
neutrino experiments, which enable a narrow parameter region to be accurately surveyed. 

In order to calculate the atmospheric neutrino intensities precisely, we need detailed 
information about (i) the primary cosmic-ray spectra at the top of the atmosphere, (ii) the 
hadronic interactions between cosmic rays and atmospheric nuclei, (iii) the propagation of 
cosmic-ray particles inside the atmosphere, and (iv) the decay of the secondary particles. 

For (i) the primary cosmic spectra, the uncertainties have been greatly reduced with new 
measurements of primary cosmic rays [7|, . Il0[|. Among these experiments, the spectra ol' 



cosmic ray protons reported by AMS and BESS show a very good agreement up to around 
100 GeV, although they were carried out in very different experimental conditions. AMS 
flew on-board the space shuttle orbiting at the altitudes between 320 km and 390 km. On 
the other hand, BESS was a balloon-borne experiment carried out at the atmospheric depths 



of about 5 g/cm^ (~ 37 km a.s.l.). Then BESS-TeV the upgraded BESS experiment 
which extended the energy region up to 540 GeV, confirmed the results of AMS and BESS. 
The results of BESS-TeV agree with that of BESS to within 3%. 

We note the observed proton spectrum by CAPRICE is obviously lower than that of AMS 
or BESS. However, the event number acquired by AMS and BESS is far larger than that 
of CAPRICE. Although it is difficult to combine the results of the different experiments, a 
combined analysis using AMS, BESS and CAPRICE data would gives a very close result 
to that using AMS and BESS only [l^. The difference of Helium spectra observed by 
AMS and BESS are also sizable. However, as the proton is the dominant component in the 
cosmic rays, the difference in terms of the nucleon flux is less than ~4% below 100 GeV. 
The nucleon flux is important in the calculation of atmospheric muon and neutrino fluxes. 
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We consider that the cosmic ray is well understood below 100 GeV, which is important for 
the calculation of atmospheric muons and neutrinos below 10 GeV. 

For the study of (ii) the hadronic interactions, the accelerator experiment is the most 
direct method. However, the data available now do not cover all the phase space necessary 
for the calculation of atmospheric neutrino flux. We study the hadronic interactions using 
the atmospheric muon flux data in this paper. As the energy of vr or K mostly goes to muons 
at their decay, the muons are considered to carry essential information of vr and K production 
in the hadronic interactions. There have been a lot of measurements of atmospheric muon 

n nn n n 

flux at ground level as compiled in Ref. pj| and at the balloon altitudes p^, |15|, |16|, 117|, 
lisl . Among them, we select the series of precise measurements of atmospheric muon s by 
BESS at various altitudes; sea level [ll|, mountain altitude 19| and balloon altitude j2o[, 
with sufficiently small systematic and statistical errors. In all these measurements, they 
used essentially the same apparatus as for the primary cosmic-ray measurements [9|, lll| . 
and systematic errors were well controlled. There are other precision measurements of the 
atmospheric muon flux useful for the study in this paper, such as the L3+C experiment j2ll |. 
Our study might be compared with the direct calculation of atmospheric neutrino flux 



from the atmospheric muon flux 



22 
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24| . in which the vr's are assumed as the dominant 



source of the atmospheric neutrinos and muons. The calculation was reviewed in Ref. 25|. 
The differences of our study from those works are in the use of a well-constructed model of 
primary cosmic ray spectra, and a full Monte Carlo simulation code for atmospheric muon 
and neutrino fluxes calculation. Then we study the hadronic interactions comparing the 
calculated and observed atmospheric muon fluxes. Both the model of primary cosmic ray 
spectra and the simulation code are the same as those used in HKKM04 calculation [l|. 
The primary flux model is very close to that constructed in Refs. 12, |2^ based on the 
AMS and BESS data. We use the primary flux model with a modification of the spectrum 
index of cosmic ray protons from -2.74 to -2.71 above 100 GeV, according to the emulsion 



chamber experiments at higher energies 27|, l28|]. The simulation code treat the (iii) propaga- 



tion of cosmic-ray particles inside the atmosphere and (iv) decay of the secondary particles 
sufficiently accurately. 

Calculating the atmospheric neutrino flux from the primary cosmic ray flux, the atmo- 



spheric muon flux was also used to calibrate the calculation in Ref. 
atmospheric density profile is crucial in this study. 



29| . However, the 
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First, we review the precision measurements of muon flux by BESS and other instruments 
in Sec. [Ill Next, we study what kind of information can be deduced from the comparison 
of the calculation of atmospheric muon flux and observed data for the calculation of atmo- 
spheric neutrino flux in Sec. IIIII We also study the effect of the atmospheric density proflle 
on the muon flux in Sec. [IV] before the comparison. Note, the seasonal change of the air 
density proflle causes ± 5 % variations of muon flux at ~ 1 GeV/c at sea level. And even 
larger variation is expected to result from changes in local meteorological conditions. In 
the same section we also discuss the effect of the uncertainty of the interaction cross-section 
between cosmic rays and air nuclei. 

We calculate the muon fluxes in the HKKM04 scheme with the observed atmospheric 
density proflle, and compared them with the precisely measured muon flux data, in Sec. El 
Note, the existence of precise atmospheric density proflle data during the observation period 
is another important reason that we use the muon flux data from the BESS measurements. 
We flnd /i"*" -|- fi~ shows reasonable agreement in the 1~30 GeV/c range between the calcu- 
lations and observations, but that the agreement worsens above 30 GeV/c. In addition, the 
u"^/u~ ratio shows a sizable difference. The difference is considered to be due to errors in 
the hadronic interaction model used in HKKM04 (DPMJET-III [30|). 

We try to improve the hadronic interaction model in Sec. IVH and compared with the data 



from recent accelerator data 



3l|. With a phenomenological consideration based on quark 



parton model, K productions are also modifled in this "improvement" . As the result of the 
modification, the observed muon fluxes are reproduced with good accuracy (Sec. IVIIj) . 

Note, the available precision muon flux data are essentially those for the vertical di- 
rections. If we have the accurately measured horizontal muon flux data, we can test the 
simulation code by the comparison of the calculated and observed muon fluxes for the hor- 
izontal directions. This comparison would be a good support for our procedure. However, 
the muon flux data for horizontal direction are poorer than those for vertical directions. We 



just show t 



rections 



le comparison of the calculation and available muon flux data for horizontal di- 



32 



331] in Sec. IVIIi The calculation of atmospheric neutrino flux with the modifled 



interaction model will be reported in the following paper 
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II. PRECISION MEASUREMENTS OF ATMOSPHERIC MUON 



The BESS group performed a series of atmospheric muon observations at various levels 
and sites; balloon altitude, mountain altitude and at sea level. 

At balloon altitudes, where the atmospheric depths (5~25 g/cm^) are much smaller than 
the interaction mean free path of protons (~ 100 g/cm^), the muon flux measurement is 
considered as an inclusive experiment with the primary cosmic ray beam and the air nu- 
cleon target. We can expect rich information about hadronic interactions from this region. 
However, only a small number of experiments had been performed, and their data are poor 
in statistics, because the muon flux itself is small at balloon altitudes and the observation 
time is limited. On January 24th, 2001, the BESS group carried out a muon flux observa- 
tion at balloon altitudes with exceptionally good statistics at Fort Sumner, NM, USA 20|. 



After reaching an altitude with a residual atmospheric depth of 5 g/cm^, the balloon slowly 
descended to 28 g/cm^ in 12.4 hours. A large number of primary and secondary cosmic rays 
were recorded during the descending period. The positive muon spectrum was obtained for 
0.50-2.55 GeV/c and negative muon spectra for 0.50-9.76 GeV/c. As discussed in Ref. 34| . 



we selected DPMJET-III as the interaction model for HKKM04 with these muon data, since 
it reproduced the observed atmospheric muon spectra better than the interaction models of 



Fritiof 1.6 



35|], Fritiof 7.02 [36|, and FLUKA'97 (37|. 



The BESS group has also measured the atmospheric muon flux for vertical directions 
at ground level; at Mt. Norikura, Japan (742 g/cm^) in September, 1999 [l9| and at 



Tsukuba, Japan (1032 g/cw?) in October, 2002 ll|. The momentum ranges covered are 
0.58-106 GeV/c at Norikura and 0.58-404 GeV/c at Tsukuba. In both experiments, the 
observation times were long enough that the systematic errors dominate the statistical er- 
rors. The overall errors were 3 % at 1 GeV/c, 3 % at 10 GeV/c, and 9 % at 100 GeV/c 
for the Norikura experiment, 2 % at 1 GeV/c, 3 % at 10 GeV/c, and 5 % at 100 GeV/c 
for the Tsukuba experiment. Note that the Tsukuba experiment was carried out with the 
BESS-TeV detector, which could not distinguish electrons and positrons from muons 11|. 
On the other hand, the BESS detector used for the Norikura experiment was equipped wit 



an electromagnetic shower counter to distinguish electrons and positrons from muons 19 1. 



As an important aspect of the muon measurement by BESS at ground level, the precise 



atmospheric density profile data are available from the Japan Meteorological Agency 
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The L3+C detector has measured the atmospheric muon flux accurately from 20 GeV/c to 
3 TeV/c and at zenith angles from 0° to 58° at CERN 2]J. The L3+C detector was originally 
constructed for the LEP experiment as the L3 muon spectrometer, and the efficiency and the 
absolute momentum scale were calibrated with the muon pairs from Z decays. The overall 
error is read from Ref. 2l| to be 4.5 % at 20 GeV/c, less than 3 % in 60 - 500 GeV/c, and 



10 % at 1.5 TeV/c. However, the L3+C detector is situated below a molasses overburden of 
30 m (6854 g/cm^), which could be a source of unknown systematic error at lower momenta. 
In the following study, we used the L3+C data for vertical directions {c^^Q zenith > 0.9), 
and in the momentum range of 60 GeV/c - 3 TeV/c. 



33| and 



For the horizontal directions, there are muon flux data from the MUTRON 
DEIS 32] experiments at sea level. MUTRON observed the muon flux from 100 GeV/c to 



20 TeV/c in momentum and from 86° to 90° (88.9° on average) in zenith angle, and DEIS 
from 10 GeV/c to 10 TeV/c in momentum and from 78° to 90° in zenith angle. However, it is 
difficult to read the systematic errors from their reports 32, |33|]. These data are potentially 
useful to study the validity of the simulation code at higher energies, rather than to study 
the hadronic interaction model. 



III. RESPONSE OF ATMOSPHERIC MUON AND NEUTRINO FLUXES TO ER- 
RORS IN HADRONIC INTERACTIONS 



Before a study of the hadronic interactions, we look at the response of atmospheric muon 
and neutrino fluxes to errors in the hadronic interaction model. If their responses are the 
same, we can study the hadronic interactions relevant to the atmospheric neutrino flux 
with the atmospheric muon flux data. Here, we derive some analytical expressions for the 
calculation of atmospheric muon and neutrino fluxes, but we actually calculate them in the 
Monte Carlo simulation, then interpret the results with the analytical expressions. 

We use the HKKM04 calculation code for the Monte Carlo calculation of atmospheric 
neutrino flux jl|. In the HKKM04 scheme, the primary cosmic ray flux model based on the 



AMS and BESS observations 



12 



26] was used with a modiflcation of the spectrum index for 



proton cosmic rays from —2.74 to —2.71 above 100 GeV, so that the extension goes through 
the center of the emulsion experiment data 271, l28| at higher energies. DPMJET-III j3J| 



was selected for the hadronic interaction model, and the US-Standard '76 



39| atmospheric 
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density profile was used. Note, however, the discussion in this section is not sensitive to the 
details of the Monte Carlo simulation scheme. 

To cut out the hadronic interactions, we write the atmospheric lepton (/z, z/^, or z/g) fluxes 

as, 

Mpi) = E / TLiPm,Pi, h) ■ Y^ip^, h) dvmdh , (1) 

where / stands for the kind of lepton, m stands for the kind of meson {tt^, K"^, ...), 
TmiPmiPh h) is the probability with which the m-meson produced with momentum and 
at altitude h creates an Z-lepton with momentum pi at ground level, and Y"^{pm,h) is the 
m-meson yield spectrum at the altitude h. As the mesons are created in the hadronic 
interaction of cosmic rays and air nuclei, the Y"^{pm, h) is written as 

Y"^{Pm, h) = pair{h) " ^ Mpproj) " VT (Pproj , Pm) " MPproj, h) dpproj , (2) 

*^ Pproj 'f^ 

where Pair{h) is the nuclear density of the air at altitude h, i stands for the kind of projectile 
{p,p, ..yH"^,..) for the hadronic interaction in the atmosphere, is the hadronic production 
cross section of the i particle and the air nuclei, f]^{Pproj,Pm) is the m-meson production 
spectrum in the hadronic interaction of i projectile and air nuclei, and (pii^Pprojjh) is the 
momentum spectrum of i particle at altitude h. Note, we assume the superposition model 
for cosmic rays heavier than protons. 

Substituting Eq. [2] in Eq. [T] and changing the integration order, we obtain 

MPi) = ^ / / / Tl:^{Pm,Pl,h)pair{h)^ai{pproj)vr{Pproj,Pm)MPproj,h)dh 

™ J Pm J Pr,rnn J H 



dpprojdpm 1 

(3) 

The term inside the square brackets in Eq. [3] is interpreted as the contribution density 
in meson production phase spaces {pproj — Pm plane) to the lepton flux. Using the scaling 
variable x = Pm./ Pproj defined in the rest frame of Air nuclei, we show the contribution density 
calculated by the Monte Carlo calculation as a scatter plot in Pproj — x plane (Fig. [1]). Note, 
the variable x is defined in the rest frame of Air nuclei. 

Here, we sampled 3,000 /i's (/i"^ + and 3,000 z/'s {y^ + t'/^ + f e + ^e) at 0.1, 1.0, 10, 
and 100 GeV/c from the HKKM04 calculation for atmospheric p and v fluxes for vertical 
directions at sea level. Then, we plotted the momenta of parent mesons, vr's (vr"'' and vr") 
or J-^'s {K"^, and Kg), and projectiles in Fig. [T]as a scatter plot. Note, here and in the 
following discussions in this section, we do not distinguish the particles from antiparticles. 
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FIG. 1: The scatter plot of vr's (vr"*" and vr", left) and ii''s {K^, and K^, right) in the phase 
space {pproj — X plane) at their production relevant to atmospheric ^u's and j^'s at the momenta, 
0.1, 1.0, 10, and 100 GeV/c, at sea level for vertical directions. Here, x is defined as 2; = Pm/Pproj 
in the rest frame of Air nuclei for m = tt or K. We sampled 3,000 /i's and 3,000 z^'s (z/g+f^) at 
each momentum from the HKKM04 calculation. The arrow in the right bottom panel shows the 
directions to which pTr,K increase. 



For the /x's, i/^'s and i/g's, most vr's or K^s are concentrated in narrow stripes for each 
momentum above 1 GeV/c. In Fig. [2l we project the points to the tt or K momentum axis 
(left) and to the x axis (right) and show them in histograms. The total number of /i's or u^s 
(z/g+z/^) is normalized to 3,000 for each momentum as in the scatter plot (Fig. [T]). However, 



the histograms for the at 100 GeV/c are multiphed by a factor of 5, due to the rapid 
decrease of z/g at this momentum. The projections to the meson momentum axis are narrow 
distributions with sharp peaks both for vr's and i^'s, which is proportional to the integrand 
of the integration in Eq. [31 The projections to x axis are very much hke each other for 
all the /i's, i/^'s, and z/g's above 1 GeV/c. The projections to x axis at 10 GeV/c are not 
shown in the figure, since they are almost the same to those at 1 GeV/c or 100 GeV/c. 
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FIG. 2: Left panel: the momentum distributions of vr's and i^'s relevant to the atmospheric /i's 
and i/'s with fixed momenta, 0.1, 1.0, 10, 100 GeV/c, at sea level for vertical directions. Right 
panel: the corresponding x distributions for all except p^^i^ = 10 GeV/c. In these figures, total 
number of /i's or z^'s (z^e+z^^) is normalized to 3,000, but the histograms are multiplied by a factor 
of 5 for Ue at 100 GeV/c. 



Changing the integration variables from dpprojdpm to [pm/ x^)dxdpmi and exchanging the 
integration order, Eq. [3] can be rewritten as. 



(W) = XI / / / Tl,{Pm,Puh)pair{h)^{a, 



[ — ,Pm)-(Pi{—^h)\dh—dpm 



X 



X 



X 



dx 



(4) 



The projection to the x axis is proportional to the integrand of the x integration, and it is 
directly connected to the hadronic interaction model. To illustrate this we introduce some 
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assumptions. First, we assume all the projectiles are nucleons. In Fig. [3l we plot the relative 
composition of the projectiles for the interactions in which the parent meson of the leptons 
are created in the HKKM04 calculation. We find the major projectiles are nucleons below 
100 GeV/c, and the contribution of meson projectile remains <15% even at higher momenta. 
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FIG. 3: The relative composition of the projectile of the last hadronic interaction, for atmospheric 
/i, Uf^, and z/e. 

Next, as the width of the Pm distributions for fixed pi is narrow in the left panel of Fig. [21 
we approximate it with the 6 function (zero width approximation) and write Tl^{pm,Pi, h) 
as, 

TmiPm,Ph h) = f^iPm, h) ■ S{pm - Pml{Ph h)), (5) 

where pm = Pmi{pi,h) is the average relation between pm, the momentum of mesons at 
altitude h, and pi, the momentum of leptons at ground. Then the pm integration is easily 
carried out and Eq. H] is rewritten as, 

MPi) = XI / / ^m(Pm, h)pair{,h)aN(^) " 07v(^, k) ■ 7]'^{'^ , Pm)dh'^dx , (6) 



X J h 



X 



X 



X^ 



with Pm = PmiiPhh). Note, as we are working with the flux sum of particles and anti- 
particles in this section, we introduced the iso-symmetric m production function: 



VN{Pproj,Pm) — {f]p {Pproj, Pm) +Vn{Pproj,Pm))/'2 



(7) 



for nucleons. Note, the argument h in the Pmi is introduced to account for the energy 
loss of fi for the leptons which are produced by fi decay. The energy loss of the mesons 
before decay is very small. The variation of the production height is estimated as ~ ± 100 
g/cm^, since the mean free path of the cosmic rays is 100 g/cm^. Therefore, the variation 
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of Pm due to the variation of production height is ~ ±0.2 GeV/c from the average, which 
is sufficiently small for leptons above 1 GeV/c. We can write the relation of pm and pi in a 
much simpler function as Pm = Pmi{pi) without any altitude dependence for pi ^ 1 GeV/c. 
Then Vi^{Pm/ ^iPrn) comes out of the h integration in Eq. [HI and it is rewritten as 



m ( Pm \ Pm 
Vn\ :Pm 



dx 



(8) 



<PiiPi) = ^ / [ / Tl^{Pm,h)pairih)aN{—)-(l)N{ — ,h)dh 

^ J X J h 

for Pi ^ 1 GeV/c. Now, the projection of contribution density to the x axis is expressed by 
the product of two terms. One stands for the hadronic interactions, and the other for the 
rest. For later convenience, we write the expression as 

(t>l{Pl) = ^<Pl{m){Pl) = / Hl^{Pm,x)dx , (9) 

m m ^ 

with Pm = Pmiipl) for / = /i, u^, and z/g and m = vr^, K"^, ■ ■ ■ . We will come back to the 
validity of the zero width approximation later. 

Here, we note it is difficult to study the hadronic interactions for K productions with 
atmospheric fi fluxes. In Fig. HI we depicted the contribution of i^'s to the atmospheric /i's 
and z/'s in the ratio to the total flux as the functions of momentum. The /^-contribution is 
limited to the atmospheric /i's below 1 TeV/c, while that the /^-contribution to atmospheric 
u is sizable above 10 GeV/c, and is dominant above 100 GeV/c for vertical directions. 
However, as most atmospheric /i's are produced in the vr-decays, we can use them to study 
the vr productions in the hadronic interactions, which are important in the calculation of 
atmospheric u flux below 100 GeV/c. 




p^^(GeV/c) 



FIG. 4: The ii'-contribution to the vertical /t's and i/'s (left), and to the horizontal /i's and z/'s 
(right). 



We continue the study of hadronic interactions concentrating on the vr productions. In 
the following, we compare the variations of the different lepton fluxes due to the error in 
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the TT production of the hadronic interaction model. This comparison should be done at 
the momentum where the parent vr's momenta are the same. Therefore, the momentum 
relation p^^ = Pniipi) is necessary for this study. We construct the function from the Monte 
Carlo data, averaging the parent vr's momenta for a fixed lepton momentum. However, 
the comparison is carried out among the fluxes of different kinds of leptons. Therefore, 
we show the momentum relation between /x's and t'^'s, and /I's and u^s as ratios, namely 
P-i^'iliPnu^iPi^^)) / P,y^ and P~^{PTTuXPye))/Pye', in Fig. O The ratio for horizontal directions is 
taken between horizontal z/^ or v^. and vertical /i. Note, the average momentum of parent 
TT for /i's are limited to > 2 GeV/c, since it is difficult for /i's produced in the decay of 
tt's with lower momenta than this limit to reach the ground level (see Fig. [2]). Therefore, 
there are no corresponding for p^^^ or p,^^ < 0.5 GeV/c. Note, the relation depends on 
tt's energy spectrum at decay, therefore, on the primary cosmic ray spectra and interaction 
model. However, changes in these do not affect the relations greatly. 
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FIG. 5: The momentum ratio of /z's to v^'s, or z/g's, whose parent vr momenta are the same on 
average. The "Horizontal" in the figure means the momentum ratio between vertical //'s and 



horizontal i/^'s or fe's- 



Let us consider that there are 2 interaction models with the hadronic vr production 
function ff^{j>-Kl x,p.„^ and rj'J^{j>.„/x,pT,), and assume they are related by the factor function 

C{p-n/x,x) as, 

(10) 



/TT/P'"' \ /-/Pt^ \ TT /Pt^ \ 

riNK — ^P^) = ■ Vn{—^P7t) 

XXX 



We will cah ri'^{pn/x,Pn) the ^-modification of ?7]^(p^/x,p,r)- The C-niodified lepton fluxes 
are calculated as, 

(11) 



[^i{n){Pi)]<; = / Ci—,x) ■ Hl{p^,x)dx , 
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where p,r = P-wiipi) for ^ = and z/g. Note, the error of the tt production function is also 

considered as a modification of the vr production function of the perfect interaction model. 

Using Eq. [TT] and the results of the Monte Carlo calculation partly shown in Fig. [21 it is 
now easy to study the effect of an error in the tt production function on the fluxes of yU, v^, 
and Vei with test functions for C,. We assume that the modification function is expanded in 



the 2nd order B-spline functions [40(] as 

a-, x) = i + Y, CU-) ■ 6(iog(x)), (12) 

i 

where the 2nd order b-spline functions use here are defined as, 

(^^ + 1.5)2 (_i.5 < < _o.5) 

^.(u) = 1 X J 1.5-2(^)2 (_o.5<^^< 0.5) 

(2^-1.5)2 (o.5<^< 1.5) (13) 

and 

Uu) =0. (^ < -1.5, or 1.5 < ^) , 

with Ui = Uq + i ■ Au for u = log(x). The B-spline function is often used to approximate a 
general continuous function. It has a compact value region and is normalized as Yli 'Ci(^) = 1- 
The differentiability is not important here. Note, when a function is approximated with the 
expansion of the B-spline functions, the variations quicker than A are suppressed. 
Then, Eq. [TT] can be rewritten as 

[0/w(pO]c = / [1 + J2^U-M^og{x))] ■ Hl{p^,x)dx , (14) 

and the difference to the original flux as, 

[S<Pi{n)iPi)]c = [<Pi{n){pi)]c - <Pi{n)iPi) = X^Q,*!— ) ■ / ^Sog{x)) ■ Hl{p^,x)dx . (15) 

j J X 

Therefore, the relative difference is calculated as, 

^hU^ = LHlip^,x)dx = 4^^^^^^-^) ■ ^^^^'^^ 

where, 

B'M = jmog{x))-Hi{p^,x)dxlj^Hl{p^,x)dx , (17) 
and p^ = P^riipi) for I = fi, z/^, and z/g. 
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The modification function corresponding to tlie error from the perfect vr production func- 
tion is expected to be a slowly varying function of x and not very different from 1 in all 
X regions. In the comparison of the vr production between an interaction model and the 
accelerator experiment data, we find typically differences of 20 ~ 30 % (see, for example. 
Fig. 15 of Ref. [261). 
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FIG. 6: The scatter plots of [5<j)u]c/4'u) for (left) and Ue (right). The variation is 

calculated for 3,000 sets of random {C^i} in [—1,1] for vertical directions at = 1 GeV/c. 

Applying the Monte Carlo data to Hi in Eq. [T71 we can consider the artificial modifica- 
tions of the TT productions function with random numbers for 100 % error in tt production 
function. Taking a set of uniform random numbers in [—1, 1] for each {C^,i}, we calculate 
the variation of [50iy^]c/0i/^ and [5</),^J ^ , using the Eq. [TBI and A = 0.5 in 

Eq. [131 The variations for 3,000 random sets are plotted as a scatter plot in Fig. 

Although the variations for the [50^] [S(pUf,]c/(pUf, or [50i/J(^/0iy<. are large, we find a 
narrow concentrations to the [S(f)fj](/(f)fj, = [S(j)iy](;/(f)iy line there. 

As C(^^i varies freely in [—1,1], the maximum difference is calculated as. 



Max 



Max <(p.)-i?r(p.) 



i 



<(p.)-i?r(p.) 

(18) 

where ly stands for z/^ and z/g. In Fig. [71 we show the maximum difference between [50/i](j/0/i 
and [S(pu^]c/4>u^, and [50/i](^/</'^ and [5(f)u^(^/ (pve as a function of p,^, and they are small for 
Pu ^ ^ GeV/c. Note, for the horizontal directions, the maximum difference between hori- 
zontal z/'s and vertical n is calculated. If we assume 20 % as the maximum error for the tt 
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production function, C^^j is sampled in [—0.2,0.2] instead of [—1,1]. Then, the maximum 
differences are multiphed by 0.2 to the values shown in Fig. [71 

For the modifications, or for the error of the hadronic tt production function, we have the 
approximate relation 

for the lepton flux whose parent momenta are the same, or Pnfi{Pij,) = Pnu^iPi^f,) = PnueiPi/J- 
Note, the approximate relation becomes invalid below ~ 1 GeV/c, as is seen from Fig. [3 

The maximum difference with A = 0.25 for the B-spline functions (Eq. [T3!) is shown with 
a dashed line in Fig. [7l However, the maximum difference with A = 0.25 is almost the same 
as that with A = 0.5, implying A = 0.5 is fine enough in this study. 
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FIG. 7: The maximum of \[S(j)uJ(/<l)u^ - and | [^^i/Jc/^i/^ - ['5</'^]c/</'^j|- The solid line 

assmnes A = 0.5 for the B-spline functions (Eq. [T3]) . and the dashed line A = 0.25. They are 
calculated analytically using the Monte Carlo data. 

We have some comments for the zero width approximation. This approximation is valid 
for the /i's until the variation of the //-energy lost in the air becomes crucial. This is because 
the energy of /i in vr-decay is limited to 



2 2 

mf. — m,, ^ - , , , „_ 



ml + ' 



2 2 

mt — mf, 

2 I 2 ' ' M ' 



mz. + m 



(20) 



or approximately in 0.73 ■ < < 1.27 ■ E^, where E^ = E.,^ ■ {m^ + m^)/2m,r is the 
average energy of /x in 7r-decay. However, the E^, in the vr-decay distributes uniformly in 
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[0, 2E^^], where Ei^^ = E,, ■ {m\ — m^)/2m^ is the average energy of E^^ in the tt decay. For 
the decay products of the /x's, the 3-body decay phase space is convoluted. Therefore, a wide 
momentum distribution is expected for the momenta of parent tt's of z/'s. Note, the /i spin 
effect is a minor effect in this discussion. However, the steep vr-decay spectrum makes the 
momentum distributions effectively narrower. Most of the momenta of parent vr's distribute 
within Ptt/S P-k ^ Pn, where is the average momentum of parent vr's. Modifying the 
delta-function in Eq. [5] to a narrow distribution function of Pm and retaining pm integration 
in the Eq. [9], we could carry out a more general study than that presented here. However, 
we expect an almost the same result due the weak dependence of the rj^{j)proj,Pm) on Pproj- 
In this section, we have studied the response of atmospheric /i and v fluxes to error in 
the TT production in the hadronic interaction model. We have shown that a modification 
affects the atmospheric /i and v fluxes originating from the tt decay at the same rate, namely 
^(pfi/'Pfj, — — '^(puj(pue: for Pu ^ ^ GeV/c. This is an important relation, since 

the error of the hadronic interaction model could be sensed by a comparison of calculated 
and observed /i flux data, especially when accurately measured fi flux data are available. 
The relation could be used not only to estimate the error in the calculation but also to tune 
the hadronic interaction model. However, this is true only when we carry out the height 
integration in Eq. [8] for /i and u correctly. In other words, the propagation of particles in 
the air must be carried out correctly (disregarding the hadronic interactions). For an error 
of the physical input which works in the same direction for the atmospheric fi and u fluxes, 
like an error in the primary cosmic ray flux model, the uncertainty may be merged in the 
uncertainty of the hadronic interaction model, and is calibrated by the atmospheric muon 
flux data collectively. However, there are some physical inputs whose error works in different 
directions for the atmospheric /i and u fluxes. We must be careful about such uncertainties. 



IV. ATMOSPHERIC DENSITY PROFILE AND INTERACTION CROSS- 
SECTION 

In this section, we study the effects of error in the atmospheric density profile and hadronic 
interaction cross section on the atmospheric muon and neutrino fluxes. Both effects are 
relatively smaller than those of hadronic interaction and primary cosmic ray flux, but the 
errors work differently on the atmospheric muon flux and neutrino flux. Therefore, it is 
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important to estimate the error for the study of hadronic interaction, then for the calculation 
of the atmospheric neutrino flux. Note, we treat the interaction cross section separately from 
the dynamics of hadronic interaction. 
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FIG. 8: Left panel: the comparison of air density profile of MSISE90[41J and US-standard '76 
in ratio at Kamioka in different seasons. The variations corresponding to e = ±5 % in Eq. [21] are 
also shown. Right panel: the same as the left panel, but at different latitude in all season average. 
The global average in the MSISE90 model is shown by the thick line. 



As the atmospheric density profile, the US-standard '76 atmosphere model [39|] is gener- 
ally used in the calculation of atmospheric neutrino flux, including the HKKM04 calculation. 
The density profile of the US-standard '76 atmosphere model is compared with that of newer 



atmosphere model MSISE90 4l| as a ratio in Fig. [HI The MSISE90 is considered the more 
realistic atmosphere model, since it gives the position and time dependent atmospheric vari- 
ations. In the left panel, we show comparison of the atmospheric density profile at Kamioka 
for different seasons. We find the maximum difference below 40 km is ~ 10 % in summer 
(Jul.) and autumn (Oct.), but that the MSISE90 air density profile is very close to that of 
US-standard '76 in winter (Jan.) and spring (Apr.). In the right panel, we show the compar- 
ison of one-year-average at different latitudes (Lat= —90, —34, 0, 34, and 90), and the global 
average with the US-standard '76. The global average agrees well with the US-standard '76 
within ~ 5 % except for very high latitude. 

To study the effect of the uncertainty of the air density profile, we consider the modifi- 
cation of the US-standard '76 air density profile as 

1 , h 



Pus,e{h) 



-■Pus[- 



(21) 
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where Pus{h) is the atmospheric density profile of the US-standard '76, and h is the altitude. 
Note the static solution for compressible gas in the gravitational field is expressed as 

p{h) = po ■ e-T^ , (22) 

and the scale height hs is proportional to the absolute temperature. Therefore, the e in 
Eq. EU corresponds to the change of the atmospheric temperature. We consider the variation 
of e = ±5 % in Eq. [2T] for the seasonal variation. Actually, the variation from winter-spring 
to summer-autumn is approximately the same as the variation of e = ~ +5 % below 
20 km (Fig. [8]). In Fig. [9], we plot the variation of atmospheric muon and neutrino fluxes as 
a ratio for the variation oi e = ±5 % in Eq. [2T1 
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FIG. 9: The variation of atmospheric muon and neutrino fluxes for a change of atmospheric density 
profile corresponding to the e ± 5 % in Eq. [21] from the US-standard '76 atmosphere model. 

Note, the variation of air density at the production height of vr's and K^s is the main 
reason for the variation of the fluxes of atmospheric neutrinos and high energy muons. The 
decay and the interaction are competitive processes for mesons (vr's and i^'s), and they 
balance when [Flight Length Before Decay] = [Interaction Mean Free Path], or 

cr ■ E/Mc" = Aair/ampNa (23) 

is satisfied. Here, r is the lifetime of the particle, E is the energy and M is the mass of 
the particle, Aair the average mass number of air nuclei, am the interaction cross section 
of the meson and air nuclei, p the mass density of the air, and Na the Avogadro constant. 
The production height is approximated by the first interaction height of cosmic rays, i.e. 
Aair/ccrNa ~ 100 g/cm^ in the air depth, where cXcr is the interaction cross section of cosmic 
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rays. In the US-standard '76 atmosphere model, it is calculated as ~ 16 km a.s.l. and the 
air density is ~ 0.16 x 10"'^ g/cm^, for vertical directions. The energies with which the decay 
and interaction balance are approximately 90 GeV for ir"^, 170 GeV for K^, and 690 GeV for 
there. Below these energies, most vr's or K^s decay producing the muons and neutrinos, 
and above these energies, most of them interact with the air nuclei, producing lower energy 
vr's or K^s. From Eq. [211 we find the air density at constant air depth decreases for 5 > 
and increases for e < 0. 

The actual production height of vr's and K^s is spread widely in the air depth, and so the 
air density there has a wide distribution. The variation of the atmospheric density profile 
changes the distribution a little, and works mildly on both neutrino and muon productions 
at higher momenta. With the variation of e = ±5 %, the neutrino flux varies ± 1.5% at 
1 GeV, ± 1.8% at 10 GeV, and ± 2.2% at 100 GeV/c. From this variation, we estimate 
that the error due to the uncertainty of atmospheric density profile would be similar to these 
values or smaller. 

The variation of muon flux at lower momenta could be explained by the change of pro- 
duction height with the muon decay. The production height approximated by the constant 
depth (~ 100 g/cm^) moves to higher altitude for e > and lower for e < 0, and so the 
muon flux increases for 5 < 0, and decreases for e > at lower momenta. The variation 
with e = ± 5%, is ~ =f5% at 1 GeV/c and larger at much lower momenta. 

Note, there are short-term variations of atmospheric density profile due to the change of 
climate corresponding to e = ± 5 % or more. They are crucial in the precise comparison of 
the calculation and observation of the atmospheric muon flux. In the following studies, we 
calculate the muon fluxes using the observed atmospheric density profile for the observation 
duration, when available. 

At the higher momenta (> 100 GeV/c), the variation of muon flux by the change of 
atmospheric density profile is smaller than the observation error of muon flux, even in the 
precision measurements. We use the US-standard '76 atmosphere model in the calculation 
for these momenta. 

A change in the cross section also changes the first interaction height, therefore, we expect 
a similar variation of lepton fluxes to the change in atmospheric density profile. However, 
the variation is a little different from that in the change of atmospheric density profile. In 
the static atmosphere model (Eq. 1221) . the mean free path of mesons at the first interaction 
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FIG. 10: The variation of atmospheric muon and neutrino fluxes resulting from a change of hadronic 
interaction cross section of cr it 5 %. 



depth is calculated as 



(24) 



in real length, where 0"^ is m-meson interaction cross section. The production depth of 
leptons are well approximated by the mean free path of cosmic ray in the column density 
calculated as Acr = Aair j o a- If we assume (Jcr/<^m ~ constant, the change of interaction 
cross section does not affect the competition of decays and interactions (Eq. [23i) . We expect 
very small variation of atmospheric neutrino flux with acr/cTm ~ constant in the calculation 
with US-standard 76 or with MSISE90. 

We study the effect of the uncertainties of the interaction cross sections ratio between 
cosmic rays and mesons, on the atmospheric muons and neutrinos. In Fig. [10] we show the 
ratios of the fluxes calculated with the variation of AcTcr = ± 5% to the flux calculated 
with the standard cross sections, keeping the am unchanged. The interaction cross section 
of nucleons are varied with that of cosmic rays. As expected, the variation of lepton flux is 
qualitatively the same with the variation due to the change of atmospheric density proflle 
shown in Fig. [9], but quantitatively, the variation is smaller than that above 1 GeV/c. The 
variation of the atmospheric neutrino flux by the change of o"cr of ±5 % is about ± 2 % at 
higher momenta, and is smaller at lower momenta. The variation of the atmospheric muon 
flux at the change of o"cr is the same at higher momenta, but is larger at lower momenta and 
in the opposite direction from that at higher momenta. 
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V. COMPARISON OF CALCULATED AND OBSERVED MUON FLUXES 
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FIG. 11: Left panel: the muon fluxes measured accurately at Tsukuba (Sept. 2002) 
Norikura (Oct. 1999) I9I by the BESS group, and at CERN by the L3-collaboration (L3+C) 
with the calculated HKKM04 muon fluxes for Tsukuba and Mt. Norikura. Right panel: the ratio 
of the muon flux data to the calculations. The dashed and dotted lines in the right panel are 
the same ratios but calculated with the US-standard '76 atmospheric model for Tsukuba and Mt. 
Norikura respectively. 



In the left panel of Fig. [HI we plot the muon fluxes measured accurately at Tsukuba 
(Sept. 2002) [ll| and on Mt. Norikura (Oct. 1999) |l9| by the BESS group, and at CERN by 
the L3-collaboration (L3-I-C) 21| with the calculated muon fluxes in the HKKM04 scheme 
for Tsukuba and Mt. Norikura. (Details of these experiments are given in Sec. [TTl) In 
the right panel of Fig. [TTl we plot the ratio [observation]/ [calculation] for those precision 
measurements to compare the calculation and observation in more detail. In the calculation 
of the muon fluxes at Tsukuba and Mt. Norikura, we used the proton and helium fluxes 
measured by the BESS group in the preceding flight carried out within 2 months to take 
into account the solar modulation of cosmic rays correctly. Also, we used the atmospheric 



density profile observed by the Japan Meteorological Agency [38| during the experimental 
periods for them. The calculations agree with the observations well (< 5%) in the range 
1~30 GeV/c. Note, the calculation with the US-standard '76 atmospheric model is also 
compared with the observed muon data in the figure. 

Below 1 GeV/c, there is a discrepancy between the two [observation]/ [calculation] ratios 
calculated for Tsukuba and Mt. Norikura. This discrepancy might be explained by a 
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different configuration of the BESS detector used for those two observations. As explained 
in Sec. Ull an electron/positron component can be distinguished from muons in the Norikura 
observation with the electromagnetic shower counter, but not in the Tsukuba observation. 
Our Monte Carlo study for the observation predicts the electron and positron production at 
the roof of the experimental hall for Tsukuba experiment, and it explain at the difference 
of calculation and observation at least qualitatively. 

In the muon observation at Norikura, surviving protons may affect the resultant muon 
flux. The [observation] /[calculatio7i\ ratio for Norikura shows some structure between 1 and 
3 GeV/c, and is systematically smaller than that for Tsukuba above 3 GeV/c. This might 
)e explained by the treatment of proton contamination in the positive muon candidates 
1^. At sea level, the proton flux is much smaller than at mountain altitude, due to the 
attenuation in the atmosphere between the two altitudes. Therefore the correction is not 
necessary for Tsukuba experiment. 

Thus, the muon fluxes in 1~30 GeV/c are well understood by the HKKM04 calculation 
with the observed atmospheric density profile. We may conclude that DPMJET-III can be 
used to calculate the atmospheric neutrino flux in 1~10 GeV region from the conclusion in 
Sec. mil However, the muon flux calculated in the HKKM04 scheme is clearly smaller than 
those observed by the precision measurements above 30 GeV/c. At those momenta, it is 
difficult to understand the difference with the uncertainty of the physical inputs, such as 
the atmospheric density profile, other than the primary cosmic ray model or the hadronic 
interaction model above 100 GeV. We note similar deficit is observed in the comparison of 
calculation with DPMJET-III and the observation of atmospheric gamma ray fiux 421 ]. 

In Fig. [121 "we show the comparison of the calculated muon charge ratio with the observed 
ones. We find the agreement of the calculation and observation are better than 10 % in the 
all energy region. However, the muon charge ratio is much more robust observation quantity, 
and refiects almost directly the tt"*"/ tt" ratio of the hadronic interaction model. The difference 
indicates an error in the charge ratio of the tt production in DPMJET-III. 

The difference of muon charge ratio at different observation levels comes from the muon 
energy loss in the atmosphere. Most of muons are produced at higher altitudes than 
Norikura, and they are observed as slightly higher momenta muons at Norikura altitudes 
than at sea level, due to the muon energy loss. 



23 



1.5 



(0 
CE 
0) 

at 



o 

c 
o 

3 



1.0 




10" 



Tsukuba (Oct. '02) 
Norikura (Sep. '99) 
L3 + C 

^ 



10 



10" 
Pli (GeV/c) 



FIG. 12: Comparison of calculation and observation for muon charge ratio. Upper data and the 
sohd line show the muon charge ratio at ground level (Tsukuba and L3+C). Lower data and the 
dashed line show the same quantities for Mt. Norikura (2770m a.s.l.). 



VI. MODIFICATION OF DPMJET-III 



Here we consider the modification of DPMJET-III [30[, without discussing the dynamics 
of the hadronic interaction, and actually apply the modification to the "inclusive DPMJET- 
III" [l|. The inclusive DPMJET-III is constructed from the output of the original DPMJET- 
III, so that it reproduces the secondary spectra of the original DPMJET-III in an inclusive 
way. In this interaction model, the conservation laws are violated in each interaction, but are 
satisfied in the statistical way. Therefore, it is not useful to reproduce an event caused by a 
single cosmic ray, but is much faster than the original interaction model. The computation 
speed is very important in the calculation of the atmospheric neutrino flux. 

Note, the region of a; > 0.1 in the it and K production is the most responsible for 
the atmospheric muons and neutrino fluxes (see Sec JIIip . We modify the gradient of the 
secondary spectra to cause changes at x ~ 0.1 for vr's and K^s without touching the multi- 
plicities. Therefore, the quantum numbers are conserved automatically. For the magnitude 
of modification, we use the ratio of the average energy before and after the modification. 

We assign a modification parameter to a valence quark of the projectile, and consider 
the same magnitude of modification for the secondary particles which have the same valence 
quark as the projectile. In p + Air interactions, the change of average energy are assigned 
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as: 



< > 


= (1 + Cu) < > 


{ud) 


< E^- > 


= (1 + q)<E0_> 


(du) 


< E^o > 
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{{uu-dd)/2) 


< Ek+ > 


= (1 + cj < > 
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< Ek- > 


= <El-> 


(sd) 


< Eko > 


= (1 + Cd) < E^o > 


ids) 


< Ej^o > 


= < El, > 


(sd) 



Here, c„ and Cd are the modification parameters assigned to the u and d quarks respectively, 
and the < Ef > is the average energy in the original DPMJET-III for the i particle. As the 
and oscillate quickly, their average energies are effectively modified as < Ej^o j^a >= 
(1 + Cd/2) < Elo >. Note, the modification of the nucleon spectra is determined after the 
modification for mesons are determined, so that the total energy is conserved to be equal to 
that of the projectile. These assumptions and parameterization naturally relate the K and 
7r productions through the parameters assigned for the u and d quarks. 

For the n + Air interactions, we assume iso-symmetry, or that the parameter for the 
(i-quark in n + Air interactions is equal to the c„ in Eq. [25l and that for the w-quark is equal 
to the Cd- As p + Air and n + Air are the major interactions in the cosmic ray propagation 
process in air, c„ and Cd are the two major parameters (Fig. [3]). 

For the energy dependence of q's, we consider polyline functions with kinks at 1, 3.16, 10, 
31.6, 100, . . . GeV. However, there are only a small number of data points above 100 GeV/c, 
and the uncertainty of the primary flux data is large there. We simply assume: 

Ci = ai- logio(£'proj/10 GeV) (26) 

above 316 GeV for z = u and d. Then, we tuned the cjs and q's at the kink points and 
Qu, ad to minimize the difference between calculations and observations. In this study, we 
used the muon flux data from L3+C at CERN above 60 GeV/c, and BESS at Tsukuba for 
all the momentum region. Note, the BESS data did not suffer from the effects of overlying 
material. 

With these procedures, we find the c„'s and q's connected by the polylines in Fig. [13] 
give the best result. The kinks seen at around 10 GeV are due to the connection to the 
NUCRIN interaction model in Figs. [T3]and[T31 Note, the differences of c„'s and c^'s result 
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FIG. 13: The best modification parameters, c^'s and c^'s, as a function of Eproj- 



from the difference in the muon charge ratio in the observation and calculation (Fig. [T2l) . 
We call thus modified interaction model as "modified DPMJET-III" . 

We compared the energy distributions to the secondary particles of the original and 
modified DPMJET-III in the left panel of Fig. [141 ^ind Z-factors in the right panel. The 
Z-factor is defined as 



and Xi 



Pi 



Pproj 



(27) 



where Ni is the multiplicity and Xi is the scaling variable used in Sec. Illll for the i secondary 
particle (7r+,7r~, ...). The power 1.7 is approximately equal to the integral spectrum index 
of the cosmic ray protons (1.71 in the primary flux model used by HKKM04). The Z-factor 
plays an essential role in the analytic calculation of the atmospheric muon and neutrino 
flux at higher energies. We flnd the Z-factors in modifled DPMJET-III have flatter energy 
dependences than the original one above 100 GeV, as is suggested by the scaling hypothesis. 

In Fig. [I5l we compared the Xp spectra of vr's production of p + Air interactions in the 
original and modifled DPMJET-III with that of p + C interactions in the NA49 experiment 
at 158 GeV/c 311] • In this comparison, we use the Feynman scaling variable deflned as 

= Pw/'^^/s in the CM-frame of a projectile and a nucleon in the target nucleus. Note, 
the scaling variable x we used in Sec. Illll is deflned a little differently from the xp using the 
momenta at the rest frame of Air nuclei. However, both deflnitions are almost equivalent 
for xp > 0.1 at 158 GeV/c. We flnd the modifled DPMJET-III reproduce the production 
spectra at x^? > 0.2 better than the original DPMJET-III. 
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FIG. 14: Left panel: The ratio of the average energy of secondary particles to the projectile energy. 
Right panel: the Z-factors for different kinds of particles. The Z-factor is defined in the text. In 
both panels, the solid fines sfiow tfie modified ones, and dashed lines show the original ones. 
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FIG. 15: The comparison of the vr's production spectra of the original and modified DPMJET-III 
in p + Air interaction with that in p + C from the NA49 experiment at 158 GeV/c. Note, the 
production spectrum of vr" of the modified DPMJET-III at 158 GeV/c is almost the same as that 
of the original DPMJETTII. 

VII. THE CALCULATIONS WITH THE MODIFIED INTERACTION MODEL 

In this section, we calculate the muon fluxes for the observations at Mt. Norikura, 
Tsukuba and CERN with the modified DPMJET-III described in the previous section. 
Except for the interaction model, we exactly repeat the calculations in Sec. |Vl i.e. we use 
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the primary flux model based on the measurement within 2 months, and observed atmo- 
spheric density proflle during the experimental period for Mt. Norikura and Tsukuba. The 
muon flux sum (//"^ + f^~) is compared in ratio [observation]/ [calculation] in the left panel 
of Fig. [ini and the muon charge ratio in the right panel. In the left panel, we also plotted 
the ratio for the horizontal muon flux data observed by the DEIS 3^ and MUTRON 33 1 
experiments. 

Comparing Fig. [16] with Figs. [TT] and [121 we flnd the agreement of calculation and ob- 
servation of the muon fluxes is greatly improved by the modiflcation in 30~300 GeV/c for 
muon fluxes at Mt. Norikura and Tsukuba and between 60 GeV~2 TeV for the CERN ex- 
periment. We summarized the remaining differences between calculations with the modified 
DPMJET-III and observations, including the experimental errors, as, 

0.04 - 0.24 ■ log(^^) p, < 1 GeV/c, 

5$^ = I 0.04 1 GeV/c <p^<20 GeV/c, and (28) 

0.04 + 0.065 ■ log(2^^^) p, > 20 GeV/c , 

and plotted these in Fig. [16] with dashed lines. We also find DEIS and MUTRON data agree 
with the calculation in the momentum ranges of 60~600 GeV/c and 200 GeV~2 TeV/c 
respectively, and are well inside the dashed lines. Note, the systematic error for DEIS and 
MUTRON are not included in the error bars. The modified DPMJET-III should be able to 
calculate the atmospheric neutrino fiux, at least for the 7r-decay, with good accuracy above 
1 GeV, from the study described in Sec. IIIII 



VIII. SUMMARY 



In this paper, we have studied the hadronic interaction for the calculation of atmospheric 
neutrino fiux, using atmospheric muon fiux data observed by precision measurements. 

We summarized the muon data from the precision measurements, and selected the data 
from BESS and L3+C for our study. There are other potentially useful data, such as muon 
observations at balloon altitudes by BESS, or for horizontal directions by MUTRON or 
DEIS. However, the former still suffers from small statistics, and the latter do not clearly 
quantify the systematic errors in their reports. 

Then, we studied the vr and K productions in the hadronic interactions of cosmic rays and 
air nuclei relevant to the atmospheric muons and neutrinos. In this study we manipulated 
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FIG. 16: The comparison of calculated muon fluxes with the modified interaction model and the 
observed ones. The dashed line is the sum of the errors in data and residuals by the modification 
(Eg. [28]) . Left panel: for the muon flux sum (/x"*" and Right panel: for charge ratio {iJ,~^ 

Note, the systematic errors for DEIS and MUTRON are not included in the error bars. 



analytic expressions, but the actual calculations were carried out in the Monte Carlo simu- 
lation using the HKKM04 calculation code. With the Monte Carlo data being interpreted 
with the use of the analytic expressions, we found the atmospheric muon and neutrino fluxes 
originated from the vr decay have the relation, A$^/<l>^ ^ A^,^^/^,^^ ~ A$i^^/$;^^ for > 1 
GeV. This relation is useful to study the error in the hadronic interaction model using the 
atmospheric muon data. 

As original DPMJET-III can reproduce the muon flux data in 1~30 GeV/c with the 
HKKM04 calculation scheme, we may say that DPMJET-III is good interaction model to 
calculation the atmospheric neutrino flux in the 1~10 GeV range. Note that vr's are the 
main source both for atmospheric muons and neutrinos in this energy region. However, the 
observed muon data show a sizable deviation from the calculated muon flux above 30 GeV/c 
and in the muon charge ratio, suggesting room for improvement in the DPMJET-III model. 
We note similar deficit is observed in the comparison of calculation with DPMJET-III and 
the observation of atmospheric gamma ray flux 421]. 

We tried an improvement of the interaction model based on the quark parton model, 
and applied the modification to the "inclusive DPMJET-III" . We have tuned the secondary 
spectra of the hadronic interactions, so that the calculation reproduces the observed muon 
fluxes accurately. As a result, the muon fluxes calculated with the modified interaction 
model agree very well with the observed muon flux data, up to ~1 TeV for both vertical 
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and horizontal directions in the flux sum (/i"*" + fi~) and up to 300 GeV/c in the charge 
ratio. The muon flux data for horizontal directions are potentially useful to examine the 
propagation code of cosmic rays in the atmosphere. However, the systematic error of the 
available muon flux data for horizontal directions ware not studied well. We just show the 
comparison of calculation in the flgure in this paper. 

The calculation of atmospheric neutrino flux, and the robustness of our modiflcation to 

n 

the DPMJET-III model, will be described in the following paper [2|]. Note, the K productions 
in the hadronic interactions are naturally modifled through the modiflcation for u and d- 
quarks. However, the modiflcation of K productions weakly couple to the muon fluxes, and 
are not tested in the comparison of calculation and observation of muon flux below ~1 TeV. 
The uncertainty of K productions and the results in atmospheric neutrino flux calculation 
is also discussed in the following paper. 
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